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We study the threshold for chaos and its relation to thermalization in the ID mean-field Bose- 
Hubbard model, which in particular describes atoms in optical lattices. We identify the threshold 
for chaos, which is finite in the thermodynamic limit, and show that it is indeed a precursor of 
thermalization. Far above the threshold, the state of the system after relaxation is governed by the 
usual laws of statistical mechanics. 
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Introduction - We study the threshold for chaos and 
the ability to thermalize of the ID mean-field Bose- 
Hubbard model (BHM) [lj. The study of thermaliza- 
tion in non-linear systems dates back to the early work 
of Fermi, Pasta, and Ulam (FPU) on a non-linear 
string, modeled by anharmonically coupled oscillators. It 
was expected that for a large number of degrees of free- 
dom, even small nonlinearities would cause the system to 
thermalize, resulting in energy equipartition. However, 
equipartition was not observed. The absence of thermal- 
ization was eventually explained in two complementary 
ways: one in terms of closeness to an integrable system, 
the Korteweg-de Vries model Q, and another in terms 
of a chaos threshold given by the theory of overlapping 
resonances put forth by Chirikov and Israilev 0, Q ■ 

Since then further studies on thermalization and ap- 
proach to equilibrium have been carried out in several 
classical field theories, including recent studies on the 
classical 4> 4 model @, Nonlinear Klein-Gordon equation 
(NLKG) 0, Non-Linear Schrodinger equation (NLSE) 
a, P , Discrete Non-Linear Schrodinger equation (DNLS) 
3, [lCj] equivalent to BHM, and Integrable Discrete Non- 
Linear Schrodinger equation (IDNLS) [i| . 

No conventional thermalization is expected in the 
NLSE and IDNLS, which are both integrable. In NLKG, 
like in FPU, the ability of the system to reach thermal 
equilibrium in the course of time evolution emerges only 
when the degree of nonlinearity exceeds a certain criti- 
cal value (see [El, [Tljj for the thermalization threshold in 
FPU) . On the contrary, the model eventually reaches 
equilibrium regardless of how small the nonlinearity is. 
In our paper we show that the BHM (along with the 
equivalent DNLS) belongs to the former class. 

Furthermore, we have compared two quantitative mea- 
sures of thermalizability: maximal Lyapunov exponent 
(whose positivity is a signature for chaos) and spectral 
entropy (which provides a distance to thermal equilib- 
rium). Both measures show a sharp threshold as one 
varies the nonlinearity strength, and the two thresholds 
are undeniably close. Furthermore, we assert that the 



chaos threshold is governed only by the parameters and 
observables that are finite in the thermodynamic limit, 
and as a result it remains finite in that limit. 

Our program is similar to a comprehensive comparison 
between FPU and </> 4 [l2j], where, however, the existence 
of the thermalization threshold in FPU is denied. 

In this paper we observe thermal behavior in time- 
averaged mean-field quantities. Note that in recent work 
on thermalization in quantum systems, thermal proper- 
ties emerge from individual quantum stationary states 
13, HI. Studies on the semi-classical regime suggest that 
the two are related, although open questions remain [l3j |. 

Empirically, our system describes the motion of 
bosonic atoms in a one-dimensional tight-binding opti- 
cal lattice [J, 

System of interest - We study the mean-field dynamics 
of an interacting one-dimensional Bose gas on a lattice 
(ID Bose-Hubbard model (BHM)) with periodic bound- 
ary conditions. The Hamiltonian in the momentum rep- 
resentation is 
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where the indices span the range n, i, j 
0. ±1, ±2, . . . , ±^f± {N s is supposed to be odd). 
Throughout the text the wavefunction ip n is normalized 
to unity: X^nl'/'"! 2 = 1- The bare frequency of each 
momentum mode is given by uo n — — 2 J cos an d 

the coupling constant is /j,q = UN a /N s . Here J and 
U are the nearest-neighbor site-hopping and on-site 
repulsion constants of the standard Bose-Hubbard 
model, respectively, and N a is the number of atoms. 
The canonical pairs are Q n — tp n , V n = iTiipni an d the 
equations of motion are given by JjVVi = — iS^F- We 
define the dimensionless non-linearity parameter, k, 
to be the ratio between the typical interaction energy 
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per site, U(N 3i /N s ) 2 , and the hopping energy per site, 
JNJN S : 
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Chaos criterion and chaos threshold from Lyapunov 
exponents - The standard signature of the chaotic nature 
of a region in phase space is that the separation between 
initially close trajectories grows exponentially with time, 
for typical trajectories, as captured by a positive max- 
imal Lyapunov exponent (MLE). In regular regions the 
separation grows linearly [16[, resulting in zero MLE. As 
we increase k in our system, we expect the phase space 
to change from being dominated by regular regions for 
small k to being dominated by chaotic regions for large 
k. In the present section, we use the MLEs to quantify 
this transition to chaos, which, as we will see in the sub- 
sequent section, coincides with a relatively broad change 
from unthermalizability to complete thermalizability. 

Consider two trajectories x(t) and x(t) with initial 
points Xo and xq, respectively. The separation 8x(t) = 
x(t) — x(t) initially satisfies a linear differential equa- 
tion, and the duration of this linear regime grows with- 
out bound as the initial separation Xq — Xq goes to zero. 
The finite-time maximal Lyapunov exponent (FTMLE) 
corresponding to the phase-space point Xq fl7i ] is 
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The limit t Rn -> oo gives the MLE, A oc (a;o). The FTM- 
LEs are themselves of intrinsic interest and in the chaotic 
regime the average over the FTMLE converges to the 
standard MLE 17, 1|| . We chose a convenient quantum 

mechanical metric, ||x — x\\ 2 = J2 n iV'n — ^J 2 ( see Gil)- 
Initially, we study the FTMLE on a 21-site lattice for 
a class of initial conditions where only the k = 0, ±1 
modes are occupied. In this subspace we sample uni- 
formly from the intersection of the microcanonical shells 
in energy and norm; the energy is chosen to be the infi- 
nite temperature energy of the subsystem, and the norm 
is 1. For each value of n, we sample 100 points, which 
we set as the initial points x . To each initial point we 
add a small random vector, as little as machine precision 
allows, to obtain the corresponding x 's. Each pair we 
propagate for a time ig n , short enough to ensure linearity 
of the evolution of Sx(t) but long enough to be able to 
clearly distinguish chaotic trajectories from regular ones 
on a plot of In Sx(t) versus t: the former increase linearly, 
and the latter, logarithmically (l8j |. We also verify that 
the average of the FTMLE 's over the ensemble of initial 
conditions does not depend on ta n as long as both criteria 
above are satisfied. In Fig. Q] the averaged FTMLEs are 
plotted as a function of the interaction strength. There is 
a distinct regime with zero Lyapunov exponent for small 
k < 0.5 and a strongly chaotic regime for k > 1 where 
all initial conditions have positive exponent. 
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FIG. 1: (color online). Averaged finite-time maximal Lya- 
punov exponent (FTMLE), A, and normalized spectral en- 
tropy, n, as functions of the nonlinearity, k. jV s =21. Inset: 
Normalized spectral entropy of final time-averaged state ver- 
sus FTMLE for each of the 100 initial condition used to com- 
pute the averaged value for k = 0.36, 0.54, 0.72, 0.9. 



Thermalizability threshold from spectral entropy - For 
coupled anharmonic oscillators, as in the FPU study, en- 
ergy equipartition among the normal momentum modes 
signified thermalization. In the BHM, the additional 
conservation of the norm modifies the quantity that is 
equipartitioned. To determine the best measure for the 
equipartition we use the variational Hartree-Fock Hamil- 
tonian [H, H HF = £„ 7ict£ p |^ n | 2 , where the set of 
Hartree-Fock energies {Twj^ f } was regarded as the varia- 
tional field. This procedure gives hcu^ F — hoj n + 2 /io Na- 
il, where /i is the chemical potential. 

The Hartree-Fock approximation is known to overesti- 
mate the interaction energy in the regime of strong inter- 
actions. For this reason, we determine the temperature 
T and the chemical potential using the time-averaged 
numerical kinetic energy (along with the norm) instead 
of the total energy. The temperature and the chemical 
potential were computed individually for each initial con- 
dition used. 

The new quantity to be equipartitioned is the 
distribution of the Hartree-Fock energy, q n (t) = 

\^ n (t)\ 2 huj^ F / J2n' |Vv0)| 2 ^' F - A quantitative mea- 
sure of the distance from thermodynamic equilibrium is 
the spectral entropy S(t) = — J2 n Qn(t) lnqjt), or more 
conveniently the normalized spectral entropy [11| , 
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where S'max = ln-ZVg is the maximum entropy, which oc- 
curs for complete equipartition of q n . In Fig. [1] the spec- 
tral entropy of the final time-averaged state, also aver- 
aged over 100 initial states (drawn from the same en- 
semble that was used for the Lyapunov exponent calcu- 
lation) is plotted for each value of k. For large nonlin- 
earities, n > 1, the normalized spectral entropy goes to 
zero, indicating remarkable agreement between the final 
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FIG. 2: (color online). Initial, final and Hartree-Fock thermal 
momentum distributions for k — 0.09, 0.45, 1.8, starting from 
the same initial state. N=21. The initial state is a represen- 
tative state and the final state is time averaged, er is the 
total energy per particle. 



state and the thermal predictions. Note that this corre- 
sponds to the chaos threshold observed previously. Fur- 
thermore, we verified that for large k, the fluctuations in 
kinetic energy scale as \^Ns, confirming their thermal na- 
ture. For k < .5 the normalized spectral entropy is above 
.5 signifying that during the time evolution the state of 
the system remains close to the initial state. As seen in 
the inset of Fig. [TJ an individual initial state with larger 
FTMLE tends to have lower spectral entropy, i.e. to relax 
to a state which is closer to the thermal one. Beginning 
at k w 0.5, where the averaged FTMLE is substantially 
non-zero, some of the initial states thermalize completely. 

In Fig. [2l the initial and time-averaged momentum 
distributions of a representative state are plotted for 
k = 0.09, 0.36 and 0.9, along with the thermal Hartree- 
Fock predictions, (\ip n \ 2 ) = (T/N a )/(Twj n + 2/j, N a - fi). 

Chaos Threshold for Different Lattice Sizes.- Let us 
start from the notion that the parameter k introduced in 
^ is the only dimcnsionlcss combination of the param- 
eters of the problem that remain finite in the thermody- 
namic limit, N B — > oo, N a /N s = const, J = const, U 
cons' . Curiously, the chaos threshold for N s = 21 is at 
k « .5, i.e. k ~ 1. Another observation comes from a re- 
lated work Q on chaos threshold in NLSE with hard-wall 
boundary conditions. The authors find that the bound- 
ary between regular and chaotic motions of momentum 
mode, n, is given by (/j,o\ijjn\ 2 ) / (nuiin) ~ 1, where hu\ is 
the lowest excitation energy, e.g. the energy of the first 
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FIG. 3: (color online). (a) Averaged Finite Time Lya- 
punov exponent, A/ J, for three different system sizes, iVs = 
11,21,41. For each k, the same energy-per-particle was used 
for each lattice size, (b) Contour lines of the Lyapunov ex- 
ponent versus the nonlinearity, k and energy-per-particle, 
£t = (H — Ho)/Na,, where H is the Hamiltonian (fTJ, and 
Ho = — 2J + (l/2)/io is the ground state value of H. N s — 11. 
The first contour line corresponds to A c = 0.02. The circles 
and dotted line give the total energies (per particle) used in 
the calculation for (a). 



excited mode in the case of the Hamiltonian (fTJ) . Assum- 
ing that the shape of the momentum distribution \ip n \ 2 
as a function of n/N s should be fixed in the thermody- 
namic limit, the left-hand-side of the above relationship 
also remains finite. These observations lead to a conjec- 
ture that the chaos criterion involves only the intensive 
parameters and observables, i.e. those that are finite in 
the thermodynamic limit. 

Our test for the above conjecture is based on the fact 
that for a chaotic motion the majority of the trajectories 
cover the whole available phase space, and as a result 
the MLE becomes, for a given set of parameters, a func- 
tion of just the conserved quantities: energy and norm. 
This implies that for the same energy-per-particle, norm, 
and nonlinearity parameter k, the Lyapunov exponents 
for different lattice sizes should be similar. In Fig. [3jt 
the averaged FTMLE is plotted for three different lat- 
tices, N s — 11, 21, and 41. For each k, the same energy- 
per-particle (in units of J) is used for all three lattices. 
The corresponding energies are shown by the solid line 
in Fig. [3)3. From the plot it is indeed evident that the 
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averaged FTMLE is universal with respect to the size of 
the lattice and that the values for N s = 11 already give 
a very good estimate of both the value of the averaged 
FTMLE and the threshold. 

Two Parametric Theory of the Chaos Threshold - The 
universality observed above suggests the most relevant 
pair of variables for mapping the chaos threshold, namely 
k and the total energy-per-particle, erf J- (In FPU one 
variable is sufficient, ultimately because there is one less 
conserved quantity.) In Fig. [3)3 contour lines of the av- 
eraged FTMLE for N s = 11 are plotted versus the non- 
linearity parameter and energy-per-particle. We use two 
sets of initial conditions with n — 0, ±1 and n = 0, ±1, ±2 
momentum modes occupied. 

One can observe a plateau in the averaged FTMLE 
for A < A c = 0.02, given by the solid line. After 
crossing the critical line the averaged FTMLE increases 
with uniform slope. The critical line resembles a hyper- 
bola with the point of closest approach to the origin at 
(k, ex) ~ (0.5, 0.2J), so that the hopping parameter J 
appears to be a relevant energy scale. This is probably 
not an accident: for ej- » J the dispersion law u n be- 
gins to deviate from the (quadratic) dispersion law of the 
integrable NLSE with periodic boundary conditions. 

Summary and outlook - In this paper we consider the 
dynamics of atoms in an optical lattice from the point of 
view of chaos theory. We identify the threshold for chaos 
and show that it corresponds to the onset of thermal- 
ization. Far above the threshold, the final state of the 
system is governed by the usual statistical mechanics. 

We see two potential applications of our results. First, 
in quantum nonequilibrium dynamics, our results can 
serve as a guide for identifying the dominant effects 
preventing thermalization in optical lattices. Based on 
the studies of the validity of the classical field the- 
ory for Bose condensates |21| our results will apply 
for the lattice site occupations satisfying N a /N s ^> 
max(«, 1) max((An/A r s ) _1 , 1), where An is the typical 
width of the momentum distribution. We note that the 
Mott regime, N a = integer x N s , An = N a , U/J > 
2.2 N a /N s [22j, lies well outside of the above criteria. 

Second, in chip-based atom interferometry with dense 
Bose condensates [23| , our results illustrate the fact that 
nonlinear instabilities cannot affect the performance of 
interferometric schemes. Recall that the force fields 
used in interferometry are usually periodic with a pe- 
riod L = A/2, where A = 2n/k, and k is the wavevector 
of light used to generate the interferometric elements. 
For spatially uniform initial conditions, the time evolu- 
tion can be described by a NLSE with periodic boundary 
conditions. In turn, the NLSE constitutes the continuum 
limit of our model, N s — > 00, where we keep constant the 
ground-state chemical potential (mq, the size of system L, 
and the ratio between the energy-per-particle Et and the 
so-called recoil energy E-r = h 2 k 2 /2m — n 2 J/N s . In this 
limit both the parameter k and the ct/J ratio tend to 



zero as N s 2 , i.e. towards the origin in Fig. [3jb), where 
the motion has no dynamical instabilities. 
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